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Higher-order gravity models have been recently the subject of much attention in the context of 
cosmic acceleration. These models are derived by adding various curvature invariants to the Einstein- 
Hilbert action. Several studies showed that these models can have late-time self-acceleration and 
could, in some cases, fit various observational constraints. In view of the infinite spectrum of invari- 
ants that could be built from curvature tensors, we propose here a method based on minimal sets of 
independent invariants as a possible route for a systematic approach to these models. We explore 
a connection made between theorems on bases in invariants theory in relativity and higher-order 
cosmological models. A dynamical system analysis is performed and some models with accelerating 
attractors are discussed. The asymptotic behavior of the models is also studied using analytical 
calculations. 

PACS numbers: 98.80.-k, 95.36.+X 



I. INTRODUCTION 

Several cosmological observations [l| established that the expansion of the universe entered a phase of acceleration. 
The cosmic acceleration can be caused by a repulsive dark energy component that has a negative equation of state, 
or, a modification to gravity at cosmological scales, see for example the reviews 0. The cosmic acceleration and the 
dark energy associated with it constitute one of the most important and challenging current problems in cosmology 
and all physics. 

As for modified gravity models, there has been recently much interest in studying higher-order gravity cosmological 
models. A bibliography search [3| shows that nearly 300 papers studying these models have appeared in the last 
3 years. The models are derived from curvature invariants that are more general than the Einstein-Hilbert action 
used in general relativity and standard cosmology, see for example 0, [H, @] and recent reviews 0, Q and references 
therein. Previous studies showed that the models can have some interesting dynamical features such as a late-time 
self-accelerating expansion without a dark energy component. It was also shown in some papers that they can have 
early-time inflation as well, e.g. [ToL [Tl|. thereby providing a unification scenario for the two cosmic accelerating 
phases. In these models, the acceleration is a consequence of how the matter is coupled to the space-time. Some of the 
models were found to fit cosmological observations [H, EH 0] and other models fit solar system tests [l^]. However, 
other papers claimed that some of the models studied fail solar system tests [16( or at least under some conditions 

El, 

but this has been contested in, for example, [18| , and also in [l9| , where models that pass these tests were discussed. 
In addition to the phenomenology, it has been argued that the models have theoretical motivations within unification 
theories of fundamental interactions and within field quantization on curved space-times 0, [24j. However, previous 
studies stressed that the models considered were chosen somewhat randomly due to the lar ge s pectrum of possible 
curvature invariants 0, [2(| and a systematic approach to these models is highly desirable [2ll |22| . Such a systematic 
approach could allow one to study methodically various possibilities in this class [U [22j and would allow one to make 
more general and more conclusive statements about the models. 

In this paper, we propose a method based on minimal sets of invariants as a possible systematic approach to the 
models. We explore an idea based on theorems from the theory of invariants in general relativity [25L 12a. l27j in order 
to develop a systematic framework that allows one to substantially reduce this very large class of model to a systematic 
list (taxonomy) of models. The idea is that curvature invariants are not inde pen dent from each other and, for a given 
algebraic type of the Ricci tensor (see for example the Segre classification [28l.l32j]) and a given Petrov type of the Weyl 
tensor (i.e. symmetry classification of space-times), e.g. |29l. l30l. l3~il. [32| . there exists a complete minimal independent 
set (basis) of these invariants in terms of which all the other invariants may be expressed. The connection between 
these theorems and higher-order cosmological models has not been explored before. As an immediate consequence of 
the proposed approach, the number of independent invariants to consider is reduced from infinity to six in the worst 
case and to only two in the cases of primary interest, including all Friedmann-Lemaitre- Robertson- Walker metrics. We 
use here the consequences of this connection in order to approach the models and study their dynamics. This should 
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set the stage for future studies where the models will be submitted to various solar and cosmological observational 
constraints. 

II. HIGHER-ORDER GRAVITY MODELS 

We refer the reader to some of the work discussing these models in some detail, see for example [1, H, @ and recent 
reviews @, Q and references therein, however, let us briefly introduce the models here. In these models, various 
curvature scalars are usually added to the Ricci scalar into the gravitational action leading to new field equations. 
The new field equations change how the content of the universe is coupled to its curvature. Let us recall that the 
Einstein field equations and the standard Friedmann cosmology are derived from variation of the Einstein- Hilbert 
action: 

at r t 

d 4 xy/^gL m (1) 

where M p ~ (8itG)(— 1/2) is the reduced Planck mass, R = g a ^R a /3 is the Ricci scalar constructed from the metric 
tensor g a ^ and the Ricci tensor R a p, and L m is the matter-energy Lagrangian. Applying variational calculus to the 
equation above leads to the usual Einstein equations and the usual Friedmann equations for standard cosmology: 

G a j3 = Ra/3 ~ 2^.9a/3 = -j^pT al3 . (2) 

The higher-order invariant models under discussion have a more general action of the form 

^f{R, R af} R a[} , R a ^ s R al37 s, R aj R a /3R^ R a(3fny R af 3jsR lS ^, ■■■) + J d 4 x^L m (3) 

where R a /3jS is the Riemann curvature tensor and / is a function of the curvature invariants. It can be quickly 
realized that there are an infinite number of curvature invariants that one can construct from combinations of the 
metric tensor, the Ricci tensor, and the Riemann tensor and their powers. 

Relevant to the question of cosmic acceleration, it has been shown in previous studies that some models with inverse 
powers of these scalars have negligible contribution to the dynamics at the proximity of strong gravitational fields 
but then start contributing significantly to the dynamics at very large (cosmological) scales producing an accelerating 
expansion. For example many papers have focused on inverse powers of R and also in some cases inverse powers of 
combinations of R, R af} R a p and R a ^ s R a01 s, e.g. 0, 

III. MINIMAL SETS OF CURVATURE INVARIANTS AND HIGHER ORDER GRAVITY MODELS 

Curvature invariants are important in general relativity studies since they allow a manifestly invariants character- 
ization of some properties of spacetimes. There has been some pioneering work by (2f| [HI, [H, HH, [H, H3, [38| and 
then a renewed interest by several recent papers [2(1 H3, Uli Hfjj. This long series of studies resulted in theorems 
about minimal sets of invariants and classifications of spacetime manifolds. It was shown for example in [26l l27j that 
there are at most 14 independent real algebraic curvature invariants in a 4-dimensional Lorentzian space. The num- 
ber of independent invariants depends on the symmetries of the spacetime as delineated by the Petrov classification 



29 , |30, |31|, |32| and also the algebraic type of the Ricci tensor as for example described by the Segre classification 
2a. |32|). The other invariants of the spacetime can be written in terms of this complete minimal set (basis). 

The Petrov classification is based on the algebraic structure of the Weyl curvature tensor C Q ^ 7 5. In an n-dimensional 
spacetime, the Weyl tensor is related to the Riemann and Ricci tensors by 

C Q( 3 7 5 = R a /3js + {gasR-fp + g^Ras - g ai Rps - g/3sR ai )/(n - 2) + R{g ai gps - g a sg/3j)/(n - l)(n - 2) (4) 

The Petrov symmetry classification types are noted by {Type I, II, D, III, N, 0} and can be understood in the 
following way. One can consider the Weyl tensor as a fourth-rank tensor operator acting on bi-vectors of space- 
time. Like in linear algebra, the eigenvalues and eigenbivectors will have some multiplicities. These multiplicities 
indicate symmetries of the Weyl tensor acting as an operator and also the symmetries of the space-time. These 
multiplicities are also related to the principal null directions and determine the Petrov type of the space-time. The 
Petrov classification's theorems give, for example, type I as the most general case with four independent null directions 
and type O as the simplest case. Current studies, using additions of invariants to the Einstein- Hilbert 's action, have 
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been done using the flat Friedmann space-times. With the approach that we propose, the simplest case (Type O) still 
covers all the Friedmann-Lemaitre-Robertson- Walker metrics. 

Next, the Segre classification is based on the algebraic structure of the trace-free Ricci tensor given by 

R 

S a /3 = Ra(3 - -£9a0, (5) 

where the tensor is viewed as a second-rank operator acting on space-time vectors. The eigenvalue equation 

S a v = Xv a (6) 

is then considered and the classification is based on the multiplicity of the eigenvalues and also whether the eigenvectors 
are null, timelike or spacelike (see e.g. Chap. 5 in [32]). For example, the Friedmann-Lemaitre-Robertson- Walker 
metrics are of Segre type Al-[(1 1 1, 1)] with two eigenvalues, one with multiplicity 1 and the other with multiplicity 
3 (three equal trace-free Ricci components) and this is the type of interest that we consider in the next section. 

It is worth clarifying that Petrov and Segre classifications depend on the symmetries of the space-time metric under 
consideration and were originally derived without reference to source fluids. Indeed, most papers as for example 
[25l . [27l . [35L l3rl l39l ] derived and discussed minimum sets of invariants in terms of the algebraic Segre and Petrov types 
without reference to source fluids, while [26] considered the type of source fluid in their discussion as obtained via 
the Einstein field equations. Since we study here modified gravity models where the field equations are not anymore 
the Einstein equations, we thus consider the minimum sets from the point of view of the algebraic Segre and Petrov 
classifications of the space-time manifold. 

Now, as it was derived in original papers on curvature invariants, the elements of the basis can be built from 
contractions of the trace-free Ricci tensor plus those of the Weyl tensor above, C a p-y6, as given by equation ((4]) and 
the complex conjugate of its self-dual tensor, i.e. Cap-^s [H, H3- 

The first element of the basis is the Ricci scalar, R. The following elements are usually put into three categories. 
The elements that are built from contraction of trace-free Ricci tensors only are noted with names starting with an r 
in the notation of 26] or with an E in the notation of as for example 

n = \E {l) =Rl = \sJS a . (7) 

Then there are pure Weyl elements starting with w in the notation of [26j or with an C in the notation of [35|, such 

as 

wi=Wl = ~C afhS C af * s , (8) 
and mixed elements starting with m in the notation of [26[ or with an D in the notation of [35[, as for example, 

mx=Ml = jC a ^ ps Sn s S aP . (9) 

We use here a notation similar to, for example, Table II in [26| except that we use Rl uppercase with no-subscript 
instead of r\ as they did. We chose this small change of notation in order to be close to the recent notation using capital 
letter for invariants used to construct higher order gravity models e.g. [|| such as P = R af3 R a /3 and Q — R a ^ 1& R a p 1 s. 

These previous studies on invariant theory were focused on building minimal sets in mathematical relativity with 
no reference to higher-order gravity models. We study here the implications of this connection with higher-order 
cosmological models and their dynamics. 



IV. HIGHER-ORDER GRAVITY MODELS BASED ON MINIMAL SETS 



Aiming to proceed in a systematic way, we consider a Petrov classification of type O (the simplest) and a Segre 
type Al-[(111),1] [H, [3^]. As mentioned earlier, this case includes all the Friedmann-Lemaitre-Robertson- Walker 
manifolds. In this case the basis of invariants reduce to only two elements, see for example [26|, |27J, |35[ , namely 

{R, Rl} (10) 

or similarly {R, ■E'(i)}- A subtle point here is that the invariant Rl = \ S a 13 Sp a is built from the trace-free Ricci tensor 
and is a true second element for the basis since it is independent from i?, unlike the commonly used P = R a ^R a fj. 
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Indeed, the trace part, R, is removed from Rl but this is not the case for P. The consequences of this are reflected 
as a significant difference between the dynamical equations that are derived respectively from models built using P 
or using Rl, in the sense that {R, Rl} provides a simpler and independent set of building blocks to construct such 
theories. In order to illustrate this, let's consider the flat Friedmann-Lemaitre-Robertson- Walker metric, 

ds 2 = -dt 2 + a(t) 2 dx 2 , (11) 

where a(t) is the scale factor. Next, we calculate the contribution of to the generalized Friedmann equation and 
find 

-2^-(H 2 -2H 2 H + 2HH) (12) 

TT4, V / V I 

where m is defined in the usual way to have dimensions of mass, an overdot signifies taking the derivative with respect 
to cosmic time, and H(t) = is the Hubble parameter used to write the third order differential equation in a(t) 

as a second order differential equation in H(t). Now, for the — *p- that was previously used in the literature, the 
addition to the Friedmann equation is given by 
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8(3H 4 + 3H 2 H + H 2 ) 3 



{H 4 + 11H 2 H 3 + 2HH 2 H + 33H 4 H 2 + 30H 6 H + 6H 3 HH + 6H S + AH 5 H). (13) 



One can see immediately that the dynamical equation from Rl is significantly reduced compared to the one built from 
P. Equations l|12p for Rl is free from any redundancy present in equation (fl3|) for P. While the two constructions 
lead to different theories, the basis {R, Rl} provides the simplest independent set of building blocks to construct 
systematically such theories. Similar simplifications are expected for other higher-order invariants models in the 
taxonomy. 

Next, we derive the general field equations for f(R, Rl) models using the action of the form 

1 = ^2 J dix ^^ R + K R > m ) + / ^x^L m (14) 
that we vary with respect to the metric and we find 

S «f> - lg>#R - ± 9 <*Pf + f R S«P + \f R g af) R + g^f R . n i - f R , a0 + 7 + \f R iS af) R 

+\{}RiS aP ).^ + \g a(3 (f R iS~<% s - i(/mS 7/ V 7 - \(fRiS ia ) ; % = SttGT^ (15) 

where we have used the definitions f R = ^ and f R i = -§^- We study the dynamics of some f(R, Rl) models in the 
next sections. 



V. MODEL is f(R,Rl) = - a ^T a ' R1)n 

Let's construct some examples based on the minimal set {R, Rl}. We use some guidance from previous studies 
[U SH, E2, S, 0, EH in order to satisfy some conditions for building models with physical degree of freedom. As 
discussed in these papers, a theory with action R + f(GB) can be re- written as the Einstcin-Hilbert action plus a 
GB-function coupled to a scalar field (f> with potential U(<p), i.e. R + f(<f>)GB — U(<f>) (45| . In the latter frame, the 
theory becomes like that of a Gauss-Bonnet one where the equations of motion (EOMs) then decouple into second 
order equations for the metric and for each scalar field involved. For example, [45| derived some conditions for GB- 
like invariants so that the EOMs do decouple into second order equations for the metric and each scalar field. It 
is worth clarifying that the decoupled equations are not the background field equations in the R + f(GB) frame 
where the Friedmann equation (the 00-field equation) contains third derivatives of the metric. As discussed in [45j |. 
the cancellation of fourth order derivatives in the decoupled equations is a necessary condition to have a theory free 
from massive spin-2 gravitons (ghosts) 0, |4l|, |4^, 0, 5|| • The condition imposes for us some relations between the 
parameters that we can use in order to combine R and Rl. The first example we consider is where the function is 
that of the pure GB invariant given by, e.g. |24l. l4ll 144 145| ■ 



GB = R 2 — AR al3 R a p + R a ^ 5 R af3lS = R 2 - AP + Q 



(16) 
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FIG. 1: LEFT: Dynamical system analysis for the case study f(R,Rl) = — r^rmzsm) ' ^ n ^ e ^ rst ort ^ er P nase portrait of 
the coordinate plane (x,v), we find an accelerating attractor with the scale factor a(t) evolving as t p where p = 5. There is 
also a repeller at p — 1. RIGHT: In the second order phase portrait with coordinates (H,H), the solid-line branches are the 
accelerating ones with attractor having the scale factor a(t) as indicated on the left. The dashed line is the separatrix solution 
where —H/H 2 = 1 and corresponds to the repeller at p = 1. The horizontal axis on the right is dH/dt. 



where we can recall the earlier conventional definitions for these scalars P — R a ^R a p and Q = R a ^ s R a pjs- Now 
using the identity 

C^C^s = C 2 = |i? 2 -2P + Q, (17) 

setting C 2 = for Friedmann manifolds, and writing the result in terms of {R, Rl} reduces the GB invariant to 

-| 8R1. We are interested in inverse powers of these invariants and, in a systematic way, the first example to 

explore should be 

f(R, Rl) = —55 . (18) 

6 

Now, to investigate if the model has a late time self-acceleration phase, we perform a dynamical system analysis. 
Following previous studies, we simplify the study to vacuum solutions to study asymptotic behavior of the universe 
since the inclusion of matter at late times as dust will not significantly contribute to the dynamics. 

However, one may question if the separatrix analysis done further below between matter dominance and vacuum 
dominance will not be affected by such a simplification. This question was addressed in [i[ where it was shown that 
adding matter to the vacuum analysis will not affect the fixed points for power-law solutions. It is also reassuring 
that some of our results on the separatrix analysis using dynamical systems for vacuum are consistent with results 
from other studies [45| that used analytical derivations and showed that such separatrix are present in similar models. 
Nevertheless, it remains of interest to verify this assumption using our approach in future work. 

Now, for the Friedmann equation from variation of R + f{R, -Rl)) reads 

3H 2 ■. (H 4 + 6H 2 H + HH + 3H 2 ) = (19) 

As usual, see for example we parameterize the second order differential equation in H(t) to a first order 
differential equation in y(x) by setting x = —H(t), y = H(t), so H = —ydy/dx by expanding the differential by 
separation of variables in dH/dt and perform a phase portrait analysis of the Friedmann equation. For cosmologically 
significant accelerating solutions, we examine power-law solutions that take the form a(t) oc t p in which we examine 
the relationship H = —H 2 /p or y = x 2 /p. In the same way, we define the asymptotic function v(x) as 

x 2 x 2 

v(x) = , ory = , (20) 

y v 
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with v 7^ so that 



dy 



-2xvdx + x 2 



(21) 



We draw attention to a strong point in if 



H 



non-zero, as v — > we anticipate the singularity 
Friedmann equation 



= \y\ — » oo, then \v\ — > if x ^ 0. Clearly, this means with x 
H — > oo. We can now convert equation (|19|) into a first order 



dv 

x— = ov 
dx 



6v z 



x 6 (72 - 2l6v + 216V 2 - 72v 3 ) 



(22) 



which allows us to study the asymptotic behavior of the solutions to this equation at late times as seen in the phase 
portrait on the left side of Fig. 1 . The attracting solutions for the asymptotic future are close to x — and at earlier 
times in the evolution of the universe lie in the negative values of x because the earlier parameterization of x = —H(t) 
denotes the x-axis as the negative of the Hubble parameter. However, since we are interested in asymptotically 
vacuum solutions, we can analytically solve for these late time solutions, as well. For this, we solve for v(x) as x — > 
in (H]) 

i) = v {5-6v Q + v 2 ) (23) 

noticing that as H — > 0, vo — » p = constant in equations (|20[) allows us to see that the asymptotic solutions of equation 
(|23jl have constant values of vq = 5,1,0 corresponding to power-law expansion with an attractor p = 5 and a repeller 
at p — 1 as seen in Fig. 1. The solution for p — vq — corresponds to y = H — > oo which is a singularity, also 
discussed in The attractor needs to be greater than 1 for accelerating expansion (this can be seen directly from 
using a(t) oc t p into the deceleration parameter ([26])), so we assume some initial conditions (a > 0) will eliminate the 
universe finding itself at solutions of p = 0, 1. In fact, these last two solutions are singular points for the generalized 
Friedmann equation (|19[) and the presence of a singularity at p = 1 is problematic, because as we discussed earlier, 
due to the equivalence of matter and vacuum fixed point solutions as shown in [4| , which we have assumed throughout 
this paper, this solution of p represents the transition point from matter domination (deceleration) to the accelerating 
phase, i.e. the point where d = 0. If this point cannot be crossed, it means that there is a separatrix at this solution, 
rendering the action cosmologically not viable. Both singularities are evident from the Friedmann equation (CG 
second order in H(t), where, the denominator equals zero for 



H + H 2 = 

where using power law solutions, a(t) = t p , translates into 

t 2 











(24) 



(25) 



and confirms the singular points p = 0,1. The singularities, if any, will come from any real solutions to equations 
equivalent to (|2"B"|) . The singular point p = 1, as confirmed by our earlier dynamical system approach was identified 
as a repeller. This is also a separatrix not allowing for the transition from matter domination to an accelerating 
phase. Where as, the attractor at p — 5 shows that the model with inverse powers of the pure GB action expressed 
in the basis {R, Rl} do have the desired accelerating late-time behavior, we should try to find actions without this 
separatrix and transition problem. 

In looking at the more traditional second-order phase portrait in the coordinate plane of (H, H) we can also see the 
attractor in Fig. 1. The solid lines leading to the accelerating attractor and the dotted line for the non-accelerating 
solution. More explicitly, in this Einstein-frame we can consider the analysis of the deceleration parameter [46j q 



-1 



H 



(26) 



where due to gravitational attraction it was thought that the universe was slowing down in its expansion at the time 
q was defined, a positive q meant deceleration. So, for q < we have acceleration as —1 — H/H 2 < 0, or we can define 
a condition for the accelerating expansion as 



(27) 
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Using this acceleration condition in the phase portrait of the right side of Fig. 1, we can see that coordinates of 
(H,H 2 ) < 1 will be leading to accelerating attractors and coordinates of (H,H 2 ) > 1 will be non-accelerating. We 
can now see in Fig. 1 that the dashed line corresponds to the separatrix, —H/H 2 = 1 because the solutions above 
this curve would be accelerating corresponding to the p = 5 attractor, and any solutions below this curve would be 
non-accelerating of which we have none for this model. Also, notice that the de Sitter solution is at (0, 1). 

In looking to eliminate other models that contain the a = separatrix, we continue in a systematic way, by 
generalizing this same analysis as actions with pure GB forms written in the basis {R, Rl} to 

f(R, Rl) = -, (28) 

where a is a dimensionless constant and n > is an integer, for both phase spaces, for which the general Fricdmann 
equations for models of this type become 

m 4n+2 

3H 2 ■. \6nH 2 H + An 2 H 2 H 

2a(2AH 2 ) n (H + H 2 ) 2 + n 

+H 2 + 2H 2 H + H A + nHH + nH 4 + n 2 HH + 2n 2 H 2 + 3nH 2 } = 0. (29) 

We can now convert equation (|29[) into a first order generalized Fricdmann equation 

dv -v 2 2 2 , v(6a — 12va + 6v 2 a)x 2 

x— = — \6nv + An 2 v - 1 + 2v - v 2 - 5n - nv 2 - An 2 ] — (30) 

alx n{n+\y m 4n+2 n24-"( z ^ +v) )~«(n + 1) 

We can analytically solve for the late-time solutions, in this general case as well. For this, we solve for v(x) as x — > 
in pO]) 

_ vq(672Vq + 4n 2 v - 1 + 2^q - vl - 5n - nv 2 - An 2 ) 

~ n(n + l) ' ' 

This asymptotic approach yields constant values of i>o = 0, 1, An + 1, where we have seen a specific study of the 
case where n = 1, (see Fig.l) allows us to determine that these solutions correspond to an accelerating attractor at 
p = An + 1, a singularity at p — and a singular repeller at p = 1. For confirmation of the separatrix, we study 
the denominator of the Friedmann equation (|29p . to find the singular points where the denominator vanishes. Again, 
even in the generalized pure GB model, the nontrivial part for this study comes from H + H 2 = which is solved in 
the same way for the singular point p = 1, as confirmed by our earlier dynamical system approach. 

So, we see that these pure GB models have a late-time accelerating power-law attractor as desired however they 
also have a separatrix and thus fail to provide a passage from matter domination to late-time self acceleration. We 
have shown this separatrix explicitly here using dynamical systems in our basis {R, -Rl}, a result that is consistent 
with previous claims (45j. In order to avoid this limitation, we will look in the next section at non-pure GB models. 



VI. MODEL 2: f(R,Rl) = 



R 2 +8R1 



Using the basis {R, Rl}, we re-express here models based on the general invariant obtained from 

aiR 2 + a 2 P + a 3 Q (32) 

where, as discussed in [45[ (see the begining of our previous section), the fourth-order derivatives in the decoupled 
equations are eliminated by using the combination a 2 — —Aa^ [45| of the dimensionless constants, oi, a 2 , and a 3 . In 
the previous section, we considered models with functions of the pure GB invariant where we also required ai = 03 but 
here we relax this condition. We use again (fT7|) with C 2 — to derive in terms of the basis {R, Rl} the combination 

ai R 2 +a 3 (--R 2 -8R1) (33) 
6 

where we can see by setting a\ = a% we return to the pure GB form. So, we must set a± ^ a 3 to avoid the separatrix 
at a — described in previous section. For example, we choose a 3 = 1 and a\ = 1/2 to build the inverse action 

mm) =w^ (34) 
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FIG. 2: LEFT: Dynamical systems for the case study f(R, Rl) = R -i/™ +8R1 ■ I n the first order phase portrait of the coordinate 
plane (x,v), we find an accelerating attractor with the scale factor a(t) evolving as t p where p ~ 2.9078 and the non-accelerating 
attractor with p w 0.96725. RIGHT: In the second order phase portrait with coordinates (H,H), the solid-line branches are 
the accelerating ones with attractor having the scale factor a(t) as indicated on the left. The dashed lines are non-accelerating. 
The horizontal axis on the right is dH/dt. 



where the minus sign in front of the function was absorbed into the denominator. Variation of this action with respect 
to the metric will lead to equations of motion that give the following Friedmann equation 

3H 2 + : (27 H 4 + 258H 2 H 3 + 744ff 4 iT 2 + 62AH 6 H 

12(3H' 2 + 8H 2 H + 8iJ 4 ) 3 

+ 128H 8 + 54:HH 2 H + 1UH 3 HH + 80H 5 H) = (35) 

The phase portrait for equation ([35]) in the coordinate space (H, H) is shown in Fig. 2, where we can see two 
attractors. With the acceleration condition we can see in the phase portrait of the right side of Fig. 2, that coordinates 
of (H, H 2 ) < 1 will be leading to accelerating attractors and coordinates of (H, H 2 ) > 1 will be non-accelerating 
attractors. The solid lines lead to the accelerating attractor, and the dashed lines lead to the non-accelerating 
attractor. Also, notice that the de Sitter solution is at (0, 1). 

For the study of power-law solutions in a(t) oc t p we convert equation (j3"5| using the earlier parameterization in 
y{x). The Friedmann equation now reads 

dv I35v 2 - 5A6v 3 + MOv 4 - 624u 5 + 128v 6 
x- 



dx 2v(27 - 72v + AOv 2 ) 

(972 - 7776w + 28512t> 2 - 59904v 3 + 76032w 4 - 55296t) 5 + 18432u 6 ; 
+ 2m 6 v(27- 72v + 40i; 2 ) 



(> 



(36) 



Its solutions in the space of (x, v) can be seen in Fig. 2. As an analytical confirmation of these results, we study 
asymptotic solutions as x — > in equation (|36|) . 



135 - 546w + 940^ - 624u 3 + 128v 4 

2(27 - 72w + 40k 2 ) ( ' 

where we eliminated the trivial solution. Solving this equation for v p = p = constant gives two real solutions and 
two complex solutions. Only the real solutions, vq = 31/16 + V 241/16 and vq — 31/16 — V241/16, are of interest 
for power-law attractors corresponding to p s» 2.9078 and p « 0.96725, respectively. Again, only the attractor p > 1 
can be of interest for late-time acceleration, so we would need to assume some initial conditions to avoid the non- 
accelerating solution. But the other attractor could serve to be useful to fit the universe to some non-accelerating 
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phases allowing for structure formation, see e.g. [4(. A check to make sure that we have avoided all the singular points 
for this model needs to be performed. A separatrix in the evolution of the universe between radiation domination, 
p = 1/2, and the transition point, p = 1 would keep the universe from reaching its accelerating phase, so we must 
make sure there are no singular points within this interval for this model. As before, we need to determine solutions 
to the denominator of the generalized Friedmann equation (|35[) 

3H 2 + 8H 2 H + 8H 4 = 0. (38) 
Now, we transform (|38|) in the same way as the previous section to write 

P^-SP + fr 2 ) = o (39) 

Solving the nontrivial part of equation (|39l) yields two complex solutions. Since there are no real non-trivial solutions, 
we have confirmed no real singular points for this model. 

Just as hoped, our results on Fig. 2 as well as the analytical confirmation show that this first model from the 
basis {i?, Rl} has a late-time accelerating phase seen as a dynamical system attractor and is free from the separatrix 
problem and thus allows for a transition from a matter dominated phase to an accelerating late-time phase. 



VII. MODEL 3: f(R, Rl) = - ( _L^_ aR1)2 

Our next example in the basis {R, Rl} is the next higher power of the invariant of last section, namely the function 

for which we derive the generalized Friedmann equation 

m 10 .... 

3H 2 : : (15H 4 + 220H 2 H 3 + 672H 4 H 2 + 5UH 6 H 

24(3iJ 2 + 8H 2 H + 8H 4 ) 4 

+64£T 8 + 60HH 2 H + 160H 3 HH + 96H 5 H) = 0. (41) 



The phase portrait for equation (|4T|) in coordinates (H, H) is plotted in Fig. 3. Once again, the condition for 
acceleration from above can be used to see which coordinate solutions (curves) will lead to the accelerating attractors 
(solid lines) and the non-accelerating solutions (dashed line). The solutions look similar to those in Fig. 2, but the 
scalings have been changed in that the solutions approach de Sitter more rapidly. The attractors (curves) have been 
pushed to the (right) H(t) axis by the more strongly driven accelerating dynamics of this model. This can be also 
seen in the following analytical examination. We analyze the Friedmann equation for these power-law solutions of 
a(t) oc t p by parameterizing it in y(x) as before. Again, the earlier definitions (|2"0|) and (|21[) are employed to express 
(ED as 



dv _ 135v 4 - 540v 5 + 864i> 6 - 544i; 7 + 6Av s 
X dx _ Av 3 (15- 40v + 24v 2 ) 

(62208u - 311040u 2 + 940032i; 3 - 1870848u 4 + 2506752w 5 - 2211840i; 6 + 1179648t> 7 - 294912w 8 - 5832)x 10 
+ 4m 10 t; 3 (15 - 40v + 2Av 2 ) 

(42) 

This first-order Friedmann equation is plotted in Fig. 3. In the phase space of (x, v) we can see our attractors for 
this system. For our analytical confirmation that these are indeed attractors for late-time acceleration, we study the 
first-order Friedmann equation (|4"2")l as x — > 0, 



135 - 5 40i;o + 864v 2 - 544^ + 64t> 4 
4(15 - 40w + 24w 2 ) 



= T " ° " TV? ° (43) 



where we again eliminate the trivial solution. This allows us to find asymptotic future attractors in v(x) phase space 
which correspond to two real solutions vq = p = constant, such that vq = 15/4 + 3\/l5/4 and vq = 15/4 — 3\/l5/4 
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FIG. 3: LEFT: Dynamical systems for the case study f(R, Rl) 



(-R'- ! /3-SRiy 2 



In the first order phase portrait of the 



coordinate plane (x,v), we find an accelerating attractor with the scale factor a(t) evolving as t p where p ~ 6.6548. There 
is also a non-accelerating solution of p ~ 0.8452. RIGHT: In the second order phase portrait with coordinates (H,H), the 
solid-line branches are the accelerating ones with attractor having the scale factor a(t) as indicated on the left. Compared to 
Fig. 2, the solid-lines here are closer to the H(t) axis nearing the de Sitter solution (0,1). The dashed lines are non-accelerating. 
The horizontal axis on the right is dH/dt. 



as well as two complex solutions. The solutions of interest for this study are power-law attractors corresponding to 
p ps 6.6548 and p « 0.8452, respectively. Again, with at least one accelerating attractor with p greater than 1, so 
we would need to assume some initial conditions to avoid the non-accelerating solution. The attractor at p ss 6.6548 
is larger than the 2.9078 obtained for the previous model confirming thus our interpretation above of the dynamical 
system plots indicating a stronger accelerating dynamics at late times for this model. We can extrapolate that even 
higher inverse powers of these first two models will have the larger of the two real solutions even more quickly approach 
the de Sitter solution. Also, upon performing a check to make sure that we do not have any singular points for this 
model we immediately notice that the nontrivial part of the denominator needed for study from the Friedmann 
equation has not changed from the previous model. The solutions that cause the denominator of the Friedmann 
equation to blow up are imaginary. Once again, we will have no real singular points for this model. Therefore, this 
model has the desired self-accelerating late-time phase and is free from the separatrix problem, allowing thus for a 
passage from matter domination to cosmic acceleration. 



VIII. MODEL 4: f(R,Rl) 



i[( ai -4)il 2 -8fll]™ 



We generalize the two previous models using the following function of the minimal set {R, Rl} and write 

4n+2 

f{R > R1) = - a [{ai -IW-8R^ (44) 

where a is a dimensionless constant, n > is an integer, and for which the general Friedmann equation reads 

m 4 "" 1 " 2 / • 

3H 2 : ■. = 1152B 2 H 6 H + 576B 2 nH s - 48BH 2 H 3 - A8BH 4 H 2 

2a6 n {6(3H 2 + 2A(BH 2 H + 24/lff 4 - H 2 ) 2 + n V 

-6nH 2 H 3 + 864(3 2 H 4 H 2 + H 4 + 288(3 2 nH 5 H + 576(3 2 n 2 H 5 H + An 2 HH 2 H + 2nHH 2 H + 792(3 2 nH 2 H 3 

+2U8f3 2 nH 4 H 2 + 2592(3 2 nH 6 H - m(3n 2 H 2 H z - 192(3n 2 H 4 H 2 - 96f3nH 2 H 3 ~ 48(3nH 3 HH + 576(3 2 H 8 

+72p 2 nHH 2 H + U4f3 2 n 2 HH 2 H - 2A(3nHH 2 H + 288(3 2 nH 3 HH - 2A[5nH 4 + 720 2 nH 4 + 36(3 2 H 4 - 12/lff 4 
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+2nH 4 - l20/3nH 4 H 2 + 2304/3 2 n 2 H 6 H + A8(3nH 5 H + 288f3 2 H 2 H 3 - 48l3n 2 HH 2 H - 96/3n 2 H 3 HH 
+5760 2 n 2 H 3 HH + 576(3 2 n 2 H 2 H 3 + lU(3nH 6 H + 2ZQ<if3 2 n 2 H 4 H 2 ) = (45) 



where we have used [3 — al — 5/6 to simplify the notation. Recalling the earlier definitions to perform the dynamical 
systems analysis, we study if actions of this general type will lead to self-acceleration at late times. The generalized 
Friedmann equation becomes 



(-36/3 2 - 8n 2 - 1 + 192/3nV - 48/3v + 288f3 2 v - 864/? V + A8/3v 2 + 72/3n - 216/3 2 n + 1 152/3 V - 6nv 



dv 
dx 

-576/3V + 96/3n 2 - 288/3V + 12/3 - 6ra - 192/3nv - 30 24/3 2 ™ 2 + 1368/3 2 nu + 2592/3W + 144/3nv 3 
+23 04/3 2 nV - 288/3n 2 u - 3456/3 2 nV - 576/3 2 nu 4 + 24f3nv 2 + 1728/3 2 n 2 v) + (& - 288/3v 2 + 288(3v 

+2 1 6/3 2 - 72/3 - 1728/3 2 , + 5184/3V - 69 1 2/3V + 3456/9V) (^) ( 6^6/3 - 24,% + 24/3, 2 - 1) 

— v 

2n(l + 48(3nv + 24/3w - 144/3 2 w + 144/3 2 w 2 + 288/3 2 nv 2 + 24/3v 2 - 288/3 2 nv - 24/3n + 36/3 2 + 72/3 2 n - 12/3 + 2n) 

(46) 

which allows us to study the asymptotic behavior of the solutions to this equation at late times as seen in the specific 
examples n = 1 and n = 2 phase portraits of Figs. 2 and 3 respectively. We are looking for late time power-law 
solutions with the form a(t) oc t p , so we need to analytically solve for the general asymptotic solutions, v(x) as x — > 
in equation (|46[) . i.e. 

= (-36/3 2 - 8n 2 - 1 + 192/3n 2 w 2 - 48/3w + 288/3 2 w - 864/3^ + 48/3t> 2 + 72/3n - 216/3 2 n + 1152/3 2 ^ - 6nv a 
-576/3 2 u 4 + 96/3ti 2 - 288/3V + 12/3 - 6n - 192(3nv a - 3024/3^ + 1368/3 2 nw + 2592/3 2 nv 3 ) + 14Apnv% 
+2304/3V?; 3 , - 288/3n 2 w - 3456/3Vv i - 576/3 2 nu 4 + 24/3nt^ + 1728/3 2 n 2 v ) X 

^va 

2n(l + A8(3nv + 2A/3v - 144/3 2 w + U4/3 2 v 2 + 288/3 2 nw 2 + 24/3v 2 - 288^ 2 tw - 24/3n + 36/3 2 + 72/3 2 n - 12/3 + 2n) 

(47) 



As usual, i? — * 0, wo — > p =constant in equations (|20|) allows us to see that the asymptotic solutions of equation (|47|) 
have the general solutions 

„ 1 V6 12/3 + 3n + 42/3n + 48^n 2 ± VS 
W0 = °'2 ± T27?' W^Tn) ' (48) 

where 3 = 240/3?i + 900/3 V + 9n 2 + 588/3n 2 + 480/3?i 3 + 2880^ 2 n 3 + 2304/3 V + 24/3. (49) 

As we will discuss below, these solutions contain at least one attractor with p > 1 that represent a cosmological model 
with a late-time self accelerating phase. Recalling that we set (3 — al — 5/6, we can see that al = 5/6 is a singular 
point to be avoided. Now, if cii — 1 we return to the pure GB form with solutions Vq = 0, 1, 4n + 1 which suffer from 
a lack of transition from matter domination phase to an acceleration phase caused by separatrix. However, we are 
interested in more general models and need to further constrain these general solutions to avoid singular points that 
could be troublesome for cosmological transitions. 

Indeed, for a model to be cosmologically viable, separatrix points (singularities) must not exist between the epoch 
for nucleosynthesis (i.e. p = 1/2) and the transition from matter domination to an accelerating phase (i.e. p = 1). 
For that, we analyze the denominator of the generalized Friedmann equation (|45[) searching for singular points. For 
power law expansion, i.e. a(t) = t p , the denominator of (|45[) vanishes for 



p 2 (6/3- 24/3p + 24/3p 2 - 1) 



t 4 

with the non-trivial solutions given by 



(50) 



^^iW (51) 
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These are also the first pair of solutions to the general dynamical system approach seen earlier in f|48|> . If the 
solutions ([51]) are real then they constitute singular points for the models. If their values are 1/2 < p < 1 then they 
represent the undesired separatrix. 

Now, (3 is negative for a\ < 5/6 and equation (|51[) has no real values thus no singular points. This will then avoid 
the separatrix problem for the generalized models of equation (|44[) . In this case, the last two solutions of equations 
(1481) represent the possible attractors. For n = 1, in the rang e 5/6 > ai > 367/507+ 10^3/169 these are negative 
attractors of no-interest. In the range 367/507+ 10 V5/169 > a x > 367/507-10^3/169 there are no real last solutions. 
And finally, in range 367/507— 10V3/169 > a\ will give two real positive attractors with one large enough for self- 
acceleration as shown in the two previous examples. In the case, n = 2, the boundaries of the three ranges change to 
5/6 > ai > 153/196 + 1^/196, 153/196 + 1^85/196 > a x > 153/196 - l-y/85/196, and 153/196 - l\/85/196 > oi, 
as separated in the n = 1 case, respectively. Finally, for the general n case, the new boundary terms can be found 
from 

ai = tt^TK? ™ T^a — 9T f 38 ™ 2 + 180 " 3 + 160n4 " 10 " - 1 ± V^f) , (52) 

3n 2 (25 + 80n + 64n 2 ) V /' y ' 

with £ = 1 + 20ri + 149n 2 + 530n 3 + 944n 4 + 800n 5 + 256n 6 , (53) 

where (+) and (— ) solutions will give ranges in the same way as cases n = 1,2. Next, for a\ > 5/6, the solutions 
Pi t 2 are real singular points and their exact values must be outside the interval [1/2, 1] in order to avoid separatrix 
problems. In this case, the previous dynamical solutions (|4"5)) are all real allowing for at least one self-accelerating 
model with p > 1. 



IX. MODEL 5: f(R, Rl) = (±i? 2 + 8R1) exp _i R2 _ SR1 



Using some guidance from previous studies of f(R) models [47j |. we continue the study of inverse functions of the 
invariant (— R 2 /3 — 8i?l) in the basis {i?, Rl} by considering the addition of an exponential function, 

f(R, Rl) = (~R 2 + 8R1) exp — , (54) 

to the action. The variation of our exponential action with respect to the metric, again yields, after some algebra, 
the Friedmann equation as 

3iJ 2 ; —. \e^(U58H 2 H 7 - 80ff 4 ij 2 + 396H i H i + 48H 6 H 3 - 2125UH 8 H i + 576H W H 

3(3,ff 2 + 8H 2 H 3 + 8H 4 ) 3 

-10368# 4 H 6 - 1008H 5 H 2 H - 1152H 7 HH - 829UH n HH - A2768H b H i H - 89856H 7 H 3 H - 432H 3 H 3 H 
-limm z H b H - U58HH 6 H - 1U048H 9 H 2 H - 9HH 2 H - 81HH 4 H - 2AH 3 HH - 83376H 6 H 5 - 64H 6 H 
+297H 2 H 5 - 27648ff 13 J ff - 16H 5 H - 2AH 2 H 3 - 300672H 10 ff 3 - 235008H 12 i? 2 - 82944i? 14 ij 
-576H 9 H + 81JJ 6 + 729H 8 + 768ff 12 )] = (55) 

for convenience we have used £ = — 6 ^ 3 f I 2 +SiI ^2f I 3 +SH 4^ ■ The phase portrait in coordinates (H,H) is not shown 
because it does not have any real solutions as we shall see by applying the phase-space parameters, in the same way 
we have used in previous models. The Friedmann equation takes the form 

x^- = -w[e s ((-440640w 4 + 466560V 5 + 21870u + 82944i; 7 - 99144w 2 - 2187 - 290304w 6 + 263088i; 3 )a; 8 
ax 

+ (2256« 5 - 1620w 4 - 576v 7 - 1152i; 6 - 81i; 2 + 567v 3 + 768v s )x i + 72v 5 + 6Av 7 - 112v 6 - 18w 4 ) + (-7128w 4 
+ 13824i; 7 - 19008w 6 - 4608w 8 - 243w 2 + 14976u 5 + 1944i; 3 )a; 6 ] / [e £: ((114048i; 4 + 27648w 6 - 11664?; + 42768t> 2 
-82944i; 5 - 89856w 3 + 1458)x 8 + (-432w 3 + 81i; 2 - 1152w 5 + 576u 6 + 1008w 4 )x 4 + 16w 6 + 9v 4 - 24w 5 )] (56) 

where £ = — - < ~ 3 ^g^ 2 +8v - ■ The phase portrait in Fig. 4 for the phase space (x,v) shows that there are no late-time 
power-law attractors. Analytically, we confirm this result with the asymptotic solution to equation (|56p as x — > 

_ 72^ + 64^ 3 -112i; 2 -18 
° _ 16v 2 + 9 - 24v (57j 
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FIG. 4: LEFT: Dynamical systems for the case study f(R, Rl) = (R 2 /3 + 8Rl) exp 
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In the first order phase portrait 



of the coordinate plane (a?,u), we see no evidence of an accelerating attractor with the scale factor a(t) evolving as t p where 
vo = p — constant. This behavior is explained in section IX. RIGHT: Dynamical systems for the case study of the expansion of 
the exponential action carried out to the first three terms f(R, Rl) = — (1 + - 



■R*/3 

the accelerating attractors from the n = 3 case dominate the series action at p — 
in section IX when n goes to infinity the solution becomes a pure de Sitter one. 



■SRI ' 2(-R^/3 

10.500 and p t 



-SRi)' 2 6(-K' ! /3-SRi) is ^ Here 
i 0.8125 however as explained 



where we eliminate the trivial solution. Which if we had power-law attractors, they would come from real solutions 
to equation (|57|) . but no real solutions exist, and thus, the exponential action does not lead to late-time power-law 
accelerating dynamics. It leads to a pure de Sitter inflationary phase. This lack of real solutions can best be understood 
by studying the behavior of the pure exponential action in a series expansion added to the original Einstein- Hilbcrt 
action. We begin with the usual series 

expz = l + z + iz 2 + iz 3 + -U 4 + 0(z 5 ) (58) 
2 b 24 

where we use z = l/(— i? 2 /3 — 8R1) and keep only up to two terms in z for simplification to derive 

f(R.Rl) = -(l+ r— — + — — !: -) (59) 

M ' V -R 2 /3-8Rl 2(-i? 2 /3-8i?l)V v ' 

We see from this series action, that at cosmological scales (long distances), the higher order term will dominate. In 
this action, the dominate term is represented by Model 4 from the previous section where n = 2, and in fact if we 
solve this action in the usual way, to find the accelerating attractors, we arrive at the two complex and two real 
solutions from Model 4, where the real solutions at vq — 15/4 ± 3\/l5/4 are the same as those plotted in Fig. 2. 
We can extrapolate here, as we have done previously, that adding higher terms to the action from the original series 
expansion, such as the third term from equation (|58p . will increase the value of n, so that the action becomes 

f(R, Rl) = - (i + _ R2 / 3 _ 8R1 + 2 (_ R 2 /3 _ 8M)2 + 6(_ jR 2/ 3 _ 8jR1 ) 3 ) ( 6 °) 

and the dominant term is represented by the n = 3 case, so that equation (|48p yields the real solutions of vq — 
181/32 ± 155/32 shown in Fig. 3. As we add more and more terms, we return to the pure exponential model above 
with an infinite number of terms, meaning n —> oo. Now, looking at the solutions given in (|48p and taking their 
limit when n — > oo, we see that in the last two solutions, the larger solution tends to infinity, or pure de Sitter 
(as for example noted for f(R) models in |47j]) and the smaller solution goes to zero. This only leaves the other 
pair of solutions there, depending only on /3, not n, which due to a\ < 5/6 give complex solutions just as we saw 
for the exponential action and the asymptotic equation (|57|) . Multiplying the infinite series by (— i? 2 /3 — 8i?l) or 
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l/(— i? 2 /3 — 8i?l) will only slow this behavior by one term, or expedite it, respectively, neither of which change the 
overall result for the exponential action. 

In sum, the analysis of the series confirms the result from Fig. 4 and equations (|57|) that the action with the 
exponential term go to a pure de Sitter behavior. 

The same techniques can be applied to study the behavior of a logarithmic action of the form (— i? 2 /3 — 
8R1) In (— i? 2 /3 — 8R1) where the expansion can be written as 



35-1 lx-1 2 lx-1 



= + +-:_^ +0 (z 4 ) (6i) 

X X o X 

which by collecting a certain number of terms can be simplified to the form 

. , s U 3 3 1 

W = T-* + **-*? (62) 

where we can see a logarithmic action of the form from (47j can be shown to have similar behavior to the exponential 
action on cosmological scales as we add more higher terms to return the series to its pure infinite form also with a 
pure de Sitter phase. 



X. CONCLUSION 



We discussed a systematic approach to higher-order gravity models based on a connection made between higher- 
order models and theorems on minimal sets of independent invariants as derived from invariants theory in general 
relativity. Various theorems demonstrate how for a given type of symmetry and source of the spacetime, there exist 
a unique basis of curvature invariant in terms of which the other invariants can be written. 

We found that the approach allows one to narrow the basic building elements to two invariants {R, Rl} in the case 
of interest covering all the Friedmann-Lemaitre- Robertson- Walker manifolds. 

We showed that as a result of the independent nature of the invariants used, the approach allows one to avoid any 
undesired redundancies in the dynamical equations. 

We investigated some examples built from the basis {R, Rl} using dynamical system analyses supplemented by 
analytical calculations and found that some combinations have cosmological evolution with late-time self-accelerating 
phase and do allow for a passage from a matter dominated phase to a cosmic acceleration phase. 

Another class of models that need to be considered using the basis approach is that of models with strong coupling 
between curvature invariants and scalar fields, see for example @, IH, [49|, [5(| and references therein. 

Finally, the models discussed here need to be subjected to further physical acceptability constraints and comparisons 
with astronomical data at solar system level and cosmological scales. Such a study will be presented elsewhere [51| . 
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